!-------------------------------------------------------------------------------------
subroutine diagnostics( dfu, Ez2D, Ex2D, f0vec, &                ! inputs 
                        t, Courant_in_z, OmegaBounce)            ! inputs
!-------------------------------------------------------------------------------------

  use dimensions

  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : kz, pi, iself, nz, runname, restart, SpeckleWidth
  use efield, only     : efield_get_epsilon, efield_get_E0z_center

  implicit none

  !-- global variables --

  real, intent(in)       :: dfu(jz, jv, jx, jDirection)
  real, intent(in)       :: Ez2D(jx,jz), Ex2D(jx,jz)
  real, intent(in)       :: f0vec(jv)
  real, intent(in)       :: t, Courant_in_z, OmegaBounce




  !-- local variables --

  logical, save          :: first_call = .true.

  real       :: PhasePhi
  real, save :: PhasePhiOld
  real       :: KE, PotE, PotEx, potEz, Phi0Energy
  real       :: nTrapped, fsq


  complex    :: epsilon
  complex    :: EzFtInternal, EzFTtotal

            
  real       :: funiform(jv)

  complex    :: Phi1Dz(jz)
  real       :: E0zCenter(jz)

  integer    :: ix, iz, iv, iDirection  ! counters


  !-----------------------------------------------------------------------------------
       
      nTrapped=0.0
      fsq=0

      call efield_get_epsilon(epsilon)        
        
      do iv=1,jv
         funiform(iv)=0.0
         do ix=1,jx
            do iz=1,jz
               do iDirection = 1,2
                  funiform(iv) = funiform(iv)+dfu(iz,iv,ix,iDirection)/(jx*jz)
                  fsq = fsq + ( f0vec(iv) + dfu(iz,iv,ix,iDirection) )**2 * dx*dz*dv
               enddo! iDir
            enddo ! z
         enddo !  x
      end do !  iv
        
        
      nTrapped = nTrapped/jz
        
      PotE=0;  PotEx=0; potEz=0; PHi0Energy=0
      do iz=1,jz
         do ix=1,jx
            PotEx = PotEx + 0.5 * Ex2D(ix,iz)**2 /(jx*jz)
            potEz = potEz + 0.5 * Ez2D(ix,iz)**2 /(jx*jz)
         enddo
      enddo
        
      potE = PotEx + potEz
        
      KE=0

      do iv=1,jv 
         KE =  KE + funiform(iv) * 0.5 * vvec(iv)**2 * dv
      enddo


      call efield_get_E0z_center(E0zCenter)

      Phi1Dz(:) = Ez2D(jx/2+1, :) - E0zCenter(:)
      call fft_forward(jz,Phi1Dz,Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFtInternal = Phi1Dz(1+nz)
        
      Phi1Dz(:) = Ez2D(jx/2 + 1, :)
      call fft_forward(jz,Phi1Dz,Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFTtotal = Phi1Dz(1+nz)

      PhasePhi = -aimag(log(EzFtInternal)) ! really phase of E
        
      do while ( ( PhasePhi - PhasePhiOld ) <  - pi )
         PhasePhi = PhasePhi + 2*pi
      enddo
      do while ( ( PhasePhi - PhasePhiOld ) >  + pi )
         PhasePhi = PhasePhi - 2*pi
      enddo
        
      PhasePhiOld = PhasePhi   

     !------------------------------
       

      if (first_call .and. restart.eq.0) then

         open(3, file=trim(runname)//'.vstime.txt')

         if(iself==1)then
            write(3, "('% 1.t  2.PhiInternal  3.PhiToTal  4.-Phase(phi)  5.KineticEnergy')")
            write(3, "('% 6.Ex^2energy  7.Ez^2energy  8.OmegaBounce  9.Courant_in_z')")
            write(3,*)
         else
            write(3, "('% 1.t   2.1/|Epsilon|  3.Re(Epsilon)  4.Im(Epsilon)  5.Re(Csi)')")
            write(3, "('% 6.Im(Csi)  7.KineticEnergy  8.Phi0Energy  9.TotalEnergy   10.f^2 ntrapped')")
            write(3,*)
         endif

      else

         open(3, file=trim(runname)//'.vstime.txt', position='append')

      end if



      if (first_call) then

         write(*,*)
         if(iself==1)then
            write(*,*)'t         Phi Internal   OmegaBounce     potEx         potEz'
         else
            write(*,*)'kt        real(Epsilon)   Total Energy     Kinetic Energy'
         endif

         first_call = .false.

      end if



      if (iself==1) then
         write(3,501) t, 2*abs(EzFtInternal) / kz, 2*abs(EzFTtotal) / kz, PhasePhi, KE, potEx, potEz, &
              OmegaBounce*SpeckleWidth*0.37233, Courant_in_z
         write(*,501) t, 2*abs(EzFtInternal) / kz, &
              OmegaBounce*SpeckleWidth*0.37233, potEx, potEz, Courant_in_z
      else
         write(3,501) t, 1/cabs(epsilon), real(epsilon), aimag(epsilon), &
              kz*kz * real( 1 - epsilon ), &
              kz*kz * aimag( 1 - epsilon ), &
              KE,  PHi0Energy, KE+PHi0Energy, fsq, nTrapped
         write(*,501) kz*t, real(epsilon), KE+PHi0Energy, KE
      endif

      close(3)

    501   format(  6400(e12.6,2x)  )


end subroutine diagnostics


!--------------------------------------------------------------------------------------
subroutine write_modulation(xMod, uPerpMod)
!--------------------------------------------------------------------------------------

  use dimensions
  use parameters, only :  runname

  implicit none

  complex, intent(in)  :: xMod(jx), uPerpMod(jx)
  integer              :: ix

  open(1, file = trim(runname)//'.modulation.txt')

  write(1,*) '% 1.x    2.xMod   3.uPerpMod'
  write(1,*)

  do ix =1, jx
     write(1,501) xvec(ix), real(xMod(ix)),  real(uPerpMod(ix))
  enddo
  close(1)


501   format(  6400(e12.6,2x)  )


end subroutine write_modulation

!--------------------------------------------------------------------------------------
subroutine write_dispersion()
!--------------------------------------------------------------------------------------

  use dimensions
  use parameters, only :  kz, runname

  implicit none

  integer              :: iz
  real                 :: phi0, phi00
  real                 :: Omega0, DeltaOmega, dOmegadPhi



  Omega0(kz) = 1.106 + 1.969*(kz-0.25) +3.16*(kz-0.25)**2 - 13.0*(kz-0.25)**3  & 
       + 12.11*(kz-0.25)**4 + 1.068*(   (0.534-kz)**0.5  -  (0.534-0.25)**0.5  )   ! accurate for 0.25<kz<0.53
 
  DeltaOmega(kz,phi0) = -(4.366 - 20.86*kz + 28.09* kz**2) * ( exp(-0.5 / kz**2) / kz**2 ) * phi0**0.5  & 
       - 343 * kz**4 * exp( -0.5 / kz**2 ) * phi0**1.5
  !  for   0.275 < kz < 0.425
 
  dOmegadPhi(kz,phi0) = -(4.366 - 20.86*kz + 28.09* kz**2) * ( exp(-0.5 / kz**2) / kz**2 ) * 0.5 / phi0**0.5  & 
       - 343 * kz**4 * exp( -0.5 / kz**2 ) * 1.5 * phi0**0.5
  !  for   0.275 < kz < 0.425

  
  
  open(1, file = trim(runname)//'.dispersion.txt')
  write(1,*)'% 1.phi   2.Vphase     3.-0.5*phixdomega/dphi'
  write(1,*)
  do iz=1,40
     phi00 = iz * 0.5/40
     write(1,501) phi00, ( Omega0(kz) + DeltaOmega(kz,phi00 ) ) /  kz, & 
          -0.5 * phi00 * dOmegadPhi(kz,phi00)
  end do
  close (1)


501   format(  6400(e12.6,2x)  )

end subroutine write_dispersion

!-------------------------------------------------------------------------------------
subroutine write_f0_initial(f0vec, dfu)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only   : runname

  implicit none

  real, intent(in)       :: dfu(jz, jv, jx, jDirection)
  real, intent(in)       :: f0vec(jv)

  real                   :: meandf
  integer                :: iv
 
  open(13, file = trim(runname)//'.f0initial.txt')
  write(13,*) '% 1.V  2.f0  3.<df>'
  write(13,*) 

  do iv=1,jv
     meandf = sum(dfu(:,iv,:,:)) / (jx*jz)
     write(13,501) vvec(iv), f0vec(iv),  meandf
  enddo

  close(13)

501   format(  6400(e12.6,2x)  )

end subroutine write_f0_initial

!-------------------------------------------------------------------------------------
subroutine write_fdist(dfu, OmegaBounce)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only : kz, vPhase, iWaveFrame, pi

  implicit none


  real, intent(in)       :: dfu(jz, jv, jx, jDirection)
  real, intent(in)       :: OmegaBounce  

 
  integer    :: ix, iz, iv
  real       :: v, f0


     f0(v)=( 1/sqrt(2*pi) )*exp(-.5*(v+vPhase*iWaveFrame)*(v+vPhase*iWaveFrame))

     open(16, file='df_vs_zv.dat')
     open(18, file='f_vs_zv.dat')
     open(13,file='f_trapped_vs_zv.dat')

     write(16,*) '% v\z';  write(16,*)  'NA', (zvec(iz),iz=1,jz)
     write(18,*) '% v\z';  write(18,*)  'NA', (zvec(iz),iz=1,jz)
     write(13,*) '% v\z';  write(13,*)  'NA', (zvec(iz),iz=1,jz)

     ix=jx/2 + 1    
     do iv=1,jv
        write(16,501) vvec(iv),  (    dfu(iz,iv,ix,1) + dfu(iz,iv,ix,2)               , iz=1,jz   ) ! dfzv
        write(18,501) vvec(iv),  (    dfu(iz,iv,ix,1) + dfu(iz,iv,ix,2) + f0(vvec(iv)), iz=1,jz   ) ! fzv
        if  (   abs( vvec(iv) ) < sqrt (2 * OmegaBounce/kz )    )then                  
           write(13,501) vvec(iv), ( dfu(iz,iv,ix,1) + dfu(iz,iv,ix,2) + f0(vvec(iv)), iz=1,jz   ) ! fzv trapped
        endif
     end do;

     close(13)
     close(18)
     close(16)

 501   format(  6400(e12.6,2x)  )


end subroutine write_fdist


!-------------------------------------------------------------------------------------
subroutine write_density_spectra(dfu, fvc0T)
!-------------------------------------------------------------------------------------

  use dimensions
  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : runname


  implicit none


  real, intent(in)       :: dfu(jz, jv, jx, jDirection)
  complex, intent(in)    :: fvc0T(jv)
 
  real       :: MeanFTvSq(jv), MeanFTzSq(jz)
  complex    :: fz(jz)
  complex    :: fvc(jv), fvc0(jv)

  integer    :: ix, iz, iv, ikz, ip
  real       :: MySum




      MeanFTvSq(:)=0
      MeanFTzSq(:)=0

      ix=jx/2
 
      do iz=1,jz 
         fvc(:) = dfu(iz,:,ix,1) + dfu(iz,:,ix,2)       
         call fft_forward(jv,fvc,fvc)
         fvc = fvc/float(jv)
         call dealias(jv,fvc)
         MeanFTvSq = MeanFTvSq + abs(fvc) **2
      end do

      do iv=1,jv         
         fz(:) = dfu(:,iv,ix,1) + dfu(:,iv,ix,2)
         call fft_forward(jz,fz,fz)  
         fz=fz/float(jz)
         call dealias(nkz,fz)
         MeanFTzSq(:)=MeanFTzSq(:)+abs(fz(:)) **2                        
      end do
 
          
      open(55,file=trim(runname)//'.DensityFTv.txt')

      write(55,"('% v-spectra in z-averaged density distribution along center line ')")
      write(55,"('% 1.p   2.<FTvSQ>   3.F0FT^2   4.Cum(p^4*<FTvSQ>) ')")
      write(55,*)
      MySum=0
      do ip=1,jv/2  - jv/6 -1
         MySum = MySum + pvec(ip)**4 * MeanFTvSq(ip)
         write(55,501) pvec(ip), MeanFTvSq(ip), abs(fvc0T(ip))**2, MySum
      end do
      close(55)
      
      open(55,file=trim(runname)//'.DensityFTz.txt')
      write(55,"('% z-spectra in v-averaged density distribution along center line')")
      write(55,"('% 1.kz   2.<FTzSQ>   3.Cum(kz^4*<FTkSQ>) ')")
      write(55,*)
      MySum=0
      do ikz=1,nkz/2  - nkz/6 -1
         MySum=MySum + kZvec(ikz)**4 * MeanFTzSq(ikz)
         write(55,501) kZvec(ikz), MeanFTzSq(ikz), MySum
      end do
      close(55)
       


    501   format(  6400(e12.6,2x)  )


end subroutine write_density_spectra


!-------------------------------------------------------------------------------------
subroutine write_umat_spectra(uzMat, uxMat)              ! inputs 
!-------------------------------------------------------------------------------------

  use dimensions
  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : u, runname, iSign

  implicit none

  real, intent(in)        :: uzMat(jx, jz, jDirection)
  real, intent(in)        :: uxMat(jx, jz, jDirection)

  integer    :: ix, iz, iv, ikz, ikx, iDirection

  real       :: avUx(jx), avUz(jx)
  complex    ::  dUx(jx),  dUz(jx)


  open(1,file=trim(runname)//".UmatFTx.txt")

  write(1, "('% x-spectra of z-averaged velocity')")
  write(1, "('% 1.kx  2.<|Ux_kx_Sq|^2>  3.<|Uz_kx_Sq|^2>  4.|dUx(kx,iz=jz/3)|^2')")
 
  do iDirection = 1, jDirection


     avUx = 0 !  real avUx(jx)    DIAGNOSTIC
     do iz=1,jz
        dUx(:) = uxMat(:,iz,iDirection)
        call fft_forward(jx,dUx,dUx)
        dUx = dUx/float(jx)
        avUx = avUx + abs(dUx)**2 / jz
     enddo  ! iz


     avUz = 0 !  real avUx(jx)    DIAGNOSTIC
     do iz=1,jz
        dUz(:) = uzMat(:,iz,iDirection)
        call fft_forward(jx,dUz,dUz)
        dUz=dUz/float(jx)
        avUz = avUz + abs(dUz)**2 / jz
     enddo

     iz = jz/3
     dUx(:) = uxMat(:,iz,iDirection) - u*iSign(iDirection)
     call fft_forward(jx,dUx,dUx)

     write(1,*)
     write(1,*)
     write(1, "('% direction  ', (I2) )" )    iDirection
     write(1,*)

     do ikx=1,nkx/2  - nkx/6 -1
        write(1,501) kxVec(ikx), avUx(ikx),  avUz(ikx),  abs(dUx(ikx))**2   
     end do

  end do ! iDirection

  close(1)



501  format(  6400(e12.6,2x)  )


end subroutine write_umat_spectra


  

!-------------------------------------------------------------------------------------
! subroutine write2D_vs_xz
!-------------------------------------------------------------------------------------

  subroutine write2D_vs_xz(fname, data)

    use dimensions

    implicit none

    real          :: data(jx,jz)
    character (*) :: fname
    

    integer ix, iz

    open(206,file=fname);

    write(206,*)'% z\x'; 
    write(206,*) 'NA', (  xvec(ix), ix=1,jx )
    do iz = 1, jz
      write(206,501) zvec(iz), ( data(ix,iz), ix=1,jx )
    enddo

    close(206)

    501   format(  6400(e12.6,2x)  )

  end

!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------


